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The thermal model of particle production is used to analyze the particle 
ratios and the p± spectra measured recently at RHIC. Our fit of the particle 
. ratios yields the temperature at the chemical freeze-out Tchem — 165 ± 7 

\^ • MeV with the corresponding baryon chemical potential /i^je„j = 41 ± 5 

I MeV. The quality of the fit shows that the model works well for RHIC. 

The p± spectra are evaluated in an approach which takes into account the 
' modifications of the initial thermal distributions by the secondary decays 

I . of resonances. All two- and three-body decays are included. This leads 

to an effective "cooling" of the spectra in the data region by about 30-40 
MeV. We find that the pion spectrum is characterized by the inverse slope 
which agrees well with the value inferred from the RHIC data. 
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1. Introduction 



Recently, a successful description of the particle ratios measured in rela- 
tivistic heavy-ion collisions has been achieved in the framework of so-called 
thermal models [1, 2, 3, 4, 5, 6, 7, 8, 9]. An important ingredient of this 
approach is the proper inclusion of the hadronic resonances, whose decays 
contribute in an essential way to the final multiplicities of the observed 
hadrons. In the first part of this paper the simple thermal model is used 
to find the optimal values of the thermodynamic parameters characteriz- 
ing the ratios measured at RHIC. In the second part we calculate the p± 
spectra of hadrons in the thermal model, including all resonances and using 
an analytic formula. The spectra, besides the original thermal distribution, 
receive large contributions from the decaying resonances. The inclusion of 
all possible decays is important, since only in that way one can obtain a 
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clean separation and identification of other effects such as expansion, flow 
[10, 11, 12, 13, 14, 15, 16], or medium modifications [17, 18, 19], which influ- 
ence the final shape of the distributions. The present analysis computes 
the spectra in a static fireball, and is introductory to the work of Ref. [20], 
where the effects of expansion are incorporated and very good agreement 
with experiment is achieved. 



2. Ratios of hadron multiplicities 

In the thermal model the particle densities are calculated from the ideal- 
gas expression 

^_9^ f°° dp 

where gi is the spin degeneracy factor of the ith hadron, Bi,Si,Ii are 
the baryon number, strangeness, and the third component of isospin, and 

Ei = ^p^ + mf. The quantities /xf^^^, /xf^em and ^^^^^^ are the chemical 
potentials enforcing the appropriate conservation laws. We note that Eq. 
(1) is used to calculate the "primordial" densities of stable hadrons and 
resonances at the chemical freeze-out. The final multiplicities receive con- 
tributions from the primordial stable hadrons, as well as from the secondary 
hadrons produced by sequential decays of resonances after the freeze-out. 
We include all light-flavor hadrons (with the appropriate branching ratios) 
listed in the newest review of particle physics [21]. We neglect the finite- 
size and excluded volume corrections.^ The temperature, Tchem^ a-nd the 
baryonic chemical potential, l^chem-: are fitted by minimizing the expression 

= ELi (^r^ - ^'■"*)^ lul , where R^'^p is the A;th measured ratio, 
(jjfc is the corresponding error, and R^^'^'^ is the same ratio as determined 
from the thermal model. The potentials Hchem a'^'^^ f-'ihem are determined 
by the two requirements: the initial strangeness of the system is zero, and 
the ratio of the baryon number to the electric charge is the same as in the 
colliding nuclei. 

Table 1 presents our fit to the particle ratios measured at RHIC. In our 
calculation, the identical ratios measured by different groups are treated 

^ If the sizes of mesons and baryons are equal, the excluded volume corrections prac- 
tically cancel in the particle ratios. This is due to the fact that the overwhelming 
majority of hadrons is heavy and may be treated as classical particles. In this case 
we may replace the Fermi-Dirac (Bose-Einstein) distribution function in (1) by the 
Boltzmann distribution function, for which the excluded volume corrections factor- 
ize. We have checked that the use of the classical statistics changes the values of the 
thermodynamic parameters by less than 3%. 
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1.02 


1.00 ± 0.02 [24], 0.99 ± 0.02[25] 




0.09 


0.08 ± 0.01 [26] 


K-/K+ 


0.92 


0.88 ± 0.05 [27], 0.78 ± 0.12 [28] 
0.91 ± 0.09 [24], 0.92 ± 0.06 [25] 


K-/TT- 


0.16 


0.15 ±0.02 [27] 


K*o/h- 


0.046 


0.060 ±0.012 [27, 29] 


iq/h- 


0.041 


0.058 ± 0.012 [27, 29] 


p/p 


0.65 


0.61 ± 0.07 [26], 0.54 ± 0.08 [28] 
0.60 ± 0.07 [24], 0.61 ± 0.06 [25] 


A/A 


0.69 


0.73 ± 0.03 [27] 




0.76 


0.82 ± 0.08 [27] 



Table 1. Our fit to the particle ratios measured at RHIC. 



separately in the definition of (number of points n = 16). In this way the 
measurements done by different groups enter independently. Very similar 
results are obtained if we first average the results of different groups to 
obtain the most likely value for each considered ratio. Our optimal value of 
Tchern = 165 ±7 MeV is consistent with the value of the critical temperature 
as inferred from the lattice simulations of QCD (with three massless flavors 
Tc = 154 ± 8 MeV, whereas with two massless flavors Tq = 173 ± 8 MeV 
[22]). We have also calculated other characteristics of the freeze-out. In 
particular, we find the energy density e = 0.5 GeV/fm^, the pressure P = 
0.08 GeV/fm^, and the baryon density pB = 0.02 fm"-^. Our calculation 
confirms the Cleymans-Redlich conjecture [5] that the energy per hadron 
at the chemical freeze-out is 1 GeV (our approach yields {E) / {N) = 1.0 
GeV). We observe, however, that the average energy per baryon is much 
larger than the average energy per meson: {Eb) / {Nb) = 1.6 GeV and 
(Em) I {Nm) = 0.9 GcV. Similar differences occur also for other heavy-ion 
collisions studied at AGS and SPS. They are caused by the different growth 
rates of the meson and baryon mass spectra [23]. In addition, we find that 
the ratios A/A and S/S are practically unaffected by the weak decays. In 
conclusion, the thermal model works well for the RHIC particle ratios, with 
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Fig. 1. The contour plot of our /n treated as a function of T^hem and M^tem- 
Our result (black circle), the result of Ref. [30] (black triangle), and the fit of Ref. 
[31] (black square) are all shown with the corresponding errors. 

thermal parameters assuming anticipated values. 

We note that our Tchem is 9 MeV lower than 174 MeV of Ref. [30], 
and 25 MeV lower than 190 MeV obtained in Ref. [31]. Nevertheless, the 
results of the three calculations are consistent within errors, as displayed in 
Fig. 1. An interesting feature shown in Fig. 1 is the characteristic valley of 
the optimal parameters. The shape of this valley indicates that the quality 
of the fit does not change if we moderately increase or decrease both Tchem 
and /xf^g^. 

3. Transverse-momentum spectra 

Next, we come to the discussion of the pT spectra. Similarly to the par- 
ticle multiplicities, the momentum spectra of the observed hadrons contain 
two contributions. The first, "primordial" contribution is purely thermal 
and is described by distribution functions characterized by the thermody- 
namic parameters at freeze-out. The second contribution comes from se- 
quential decays of the resonances. Since a substantial part of the produced 
particles comes from the decays, one may expect that the measured spectra 
are significantly changed by this effect. Indeed, it has been already known 
for a long time that the resonance decays modify the \ow-pT spectrum, due 
to the limited phase space of the emitted particles [32, 33, 34]. So far, how- 
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ever, a rather hmited number of the resonances has been included in such 
analyses. We take into account all two- and three-body hadronic decays. 
In this study we neglect other effects which can change the spectra, such as 
expansion or flow. 

Below we sketch our method of dealing with two- and three-body decays. 
The initial momentum distribution function of a resonance in the fireball rest 
frame, f{k), as well as the momentum distribution of the emitted particle in 
the resonance's rest frame, are isotropic. The latter is true, since we average 
over all possible polarization states. Thus, in the case of a two-body decay, 
the spectrum of the emitted particle (denoted by index 1) is obtained from 
the expression 

where p is the momentum of the emitted particle in the rest frame of the 
resonance, and is the Lorentz transformation to the fireball rest frame, 

LkV = V> + [{lk-'^) vl-Vk-V + lk E*\vk- (3) 



Here k is the momentum of the resonance in the fireball, = k/y A;^ + m^, 
and 7fe = (1 — f|) '^^'^ . In Eq. (2) p* is the magnitude of p and E* = 



p1 + m\ . The standard formula gives 

_ - (mi - m2f){m\ - (mi + ma)^)] ^'^ , . 

where m/j is the mass of the resonance, whereas mi and m2 are the masses 
of the emitted particles. The quantity h is the branching ratio for the con- 
sidered channel, and Jr and Ji are the spins of the resonance and particle 
1, respectively. The physical interpretation of Eq. (2) is clear: the isotropic 
distribution of particle 1 in the resonance rest frame, 5(|p| — p*) /(47rp*^), 
is boosted to the fireball frame, and there folded with the resonance distri- 
bution f{k). 

In order to write Eq. (2) in a more compact form, we do the change of 
variables: 

p' = L,p, dh = ^dy = ^dy. (5) 

The integration over p' in (2) becomes trivial and we find 
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Since the quantity L^^<1 is the momentum of the emitted particle in the 
reference frame connected with the resonance, we may rewrite Eq. (6) in 
the explicitly Lorentz-covariant way 

/i(|q|) = ^ !t^! T^^ [d'kf{k)5('^-EA. (7) 
•"^'^'^ 2Ji + 1 4.Trp* Eg J ^ ^ \ ruR J 

Here k'^ and are the four-momenta of the resonance and of the emit- 
ted particle, respectively. We note that Eq. (7) was used previously by 
Sollfrank, Koch, and Heinz [32]. In their approach, however, no assump- 
tions about the isotropic emission were made, so f{k) and consequently fi{q) 
were treated as the functions of two arguments: rapidity and transverse mo- 
mentum. In our approach f(k) and fi{q) depend only on the magnitude of 
the three-momenta, and integration over the polar and azimuthal angles in 

(7) can be done analytically. This leads to a simple expression 

f'^'^^ = ' 2J, + l2Egp*q L,,, '''f^'^^ 

where the limits of the integration are k±{q) = mn \E*q ^p*Eq\ /m^. Eq. 

(8) is a relativistic generalization of the formula derived in Ref. [34]. Its 
simplicity turns out to be especially important in the numerical treatment 
of the sequential decays. 

In the case of three-body decays we can follow the same steps as above, 
with the extra modification connected with the fact that different values 
of p* are possible now. This introduces an additional integration in Eq. 
(2). The distribution of the allowed values of p* may be obtained from the 
phase-space integral 

J -C/pl -C^p2 

where Pi,P2 and p3 arc the momenta of the emitted particles, Ep-^.Ep.^ 
and Ep^ are the corresponding energies (all measured in the resonance rest 
frame), M is the matrix element describing the three-body decay, and N 
is the normalization constant. For sake of simplicity we assume, similarly 
as in [32], that M can be approximated by a constant. Operationally, the 
final expression for three-body decays is a folding of two-body decays over 
p* with a weight following from elementary considerations based on Eq. (9) . 

Our numerical procedure is as follows: We initialize the spectra of all 
hadrons as given by the thermal model. Next, we start from the heaviest 
particle, and proceed with the decay, thus feeding the spectra of the prod- 
ucts. We repeat this step for all particles, going down with the mass, untill 
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Fig. 2. The transverse-mass distributions of pions, kaons, protons and antiprotons. 

all resonances are taken into account. Our results for the p± spectra are 
shown in Fig. 2. The initial thermal distributions are represented by the 
dashed lines, whereas the final distributions are represented by the solid 
lines. In the case of pions the long-dashed line shows the result of including 
two-body decays only. We observe that the final distributions of hadrons 
are considerably steeper than the original distributions. We term this phe- 
nomenon as "cooling" of the spectra by secondary decays. For the pions the 
initial inverse slope, calculated in the range 0.3 GeV < p± < 0.9 GeV [35], 
equals 219 MeV. Hadronic decays cool it by 34 MeV down to 185 MeV. The 
last value agrees well with the measurement of the PHENIX Collaboration 
[35] . ^ The inverse slopes obtained for other hadrons are smaller than those 
inferred from the data. Definitely, other processes influence the observed 
spectra. 

^ Contrary to a naive expectation, tlie inverse slopes of the thermal distributions in 
Fig. 2 (dotted lines) do not correspond to the physical temperature characteriz- 
ing the thermal distribution /. This effect is induced by the presence of m± as 
a prefactor multiplying the thermal distribution in the expression for the yield: 
dN / (mTdmTdy) ^ m±cosh{y)f{m±cosh{y)). Since the values of m± are not asymp- 
totic, the inverse slopes obtained in the region of p± ~ 0.5 — 1 GeV are considerably 
higher than the temperature. 
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Our calculation does not show any peaked low-p± enhancement in the 
pion spectrum. It would, if for instance we had included only the A(1232) 
decays. However, with all decays included, the increase of the spectrum is 
uniform due to the fact that different decays populate different momenta. 
The shape of the pion spectrum in Fig. 2 is concave, thus the thermal 
model does not reproduce the experimental convex shape [35] . One can see 
that the contributions from the thrcc-body decays to the pion spectrum are 
important at small momenta, however, they are smaller than the contribu- 
tions from the two-body decays. The impact of three-body decays on other 
spectra is negligible. 

Certainly, the expansion or flow effects further modify the spectra. Their 
inclusion is necessary to obtain the complete agreement with the data [20]. 
However, the effect of "cooling" from secondary decays described in this 
paper is an important ingredient of any analysis of momentum spectra in 
relativistic heavy-ion collisions. 
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